#include "mcmcFaceModel.h"


#define COMPLEXTY_FACE		0.10


// operators
MCMC_FaceRect & MCMC_FaceRect::operator=(const MCMC_FaceRect &a)
{
	if (this != &a)
	{
		rect_face	= a.rect_face;
		energy		= a.energy;
		vc_a_face	= a.vc_a_face;
		dmin		= a.dmin;
		dmax		= a.dmax;
	}

	return *this;
}


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

MCMC_FacePCA::MCMC_FacePCA()
{
}

MCMC_FacePCA::~MCMC_FacePCA()
{
}

/* PCA FACE
       
	face1   face2  face3
	1		1		1
	2		2		2
	.       .       .
	nxn     nxn    nxn
*/
bool MCMC_FacePCA::LoadFromFile(const char *str_name)
{
	ifstream file;

	file.open(str_name, ios::in|ios::binary|ios::_Nocreate);
	if (file.fail())
		return false;


	int	size,num_component;

	file.read((char*) &size, sizeof(size));
	file.read((char*) &num_component, sizeof(num_component));

	int j;

	mx_PCAs.SetDimension(size, num_component);
	for (j=0; j<mx_PCAs.rows(); j++)
		file.read((char*)&(mx_PCAs(j, 0)), sizeof(mx_PCAs(j, 0))*mx_PCAs.cols());

	vc_sigmas.SetRowVector(num_component);
	file.read((char*)&(vc_sigmas(0)), vc_sigmas.size()*sizeof(vc_sigmas(0)));
	
	vc_mean.SetRowVector(size);
	file.read((char*) &vc_mean(0), vc_mean.size()*sizeof(vc_mean(0)));

	file.read((char*) &de, sizeof(de));

	file.read((char*)&scale_m, sizeof(double));

	file.read((char*)&scale_s, sizeof(double));

	file.close();

	return true;
}


void MCMC_FacePCA::Projection(mcmcVector<double> &vc_a, const mcmcVector<double> &vc_V)
{
	assert(vc_V.size() == SizeOfPCA());

	vc_a.SetRowVector(NumComponents());

	int i,j;

	for (i=0; i<vc_a.size(); i++)
	{
		vc_a(i) = 0.0;

		for (j=0; j<vc_V.size(); j++)
			vc_a(i) += mx_PCAs(j, i)*(vc_V.Data(j)-vc_mean(j));
	}
	

	mcmcVector<double> vc_temp;
	vc_temp.SetRowVector(vc_V.size());
	vc_temp.InitValue(0.0);

	for (i=0; i<NumComponents(); i++)
		for (j=0; j<SizeOfPCA(); j++)
			vc_temp(j) += vc_a(i)*mx_PCAs(j, i);

	for (j=0; j<SizeOfPCA(); j++)
		vc_temp(j) += vc_mean(j);
}

void MCMC_FacePCA::Reconstruct(mcmcVector<double> &vc_V, const mcmcVector<double> &vc_a)
{
	mcmcVector<double> vc_temp;

	vc_temp = vc_a;

	for (int i=0; i<vc_a.size(); i++)
		vc_temp(i) *= (vc_sigmas(i) / (vc_sigmas(i) + de));

	vc_V.SetRowVector(SizeOfPCA());
	vc_V.InitValue(0.0);

	for (int j=0; j<SizeOfPCA(); j++)
	{
		for (int i=0; i<NumComponents(); i++)
			vc_V(j) += mx_PCAs(j, i)*vc_temp(i);
		vc_V(j) += vc_mean(j);
	}
}

void MCMC_FacePCA::Reconstruct(Raster<double> &ras_image, const mcmcVector<double> &vc_a, double dmin, double dmax)
{
	mcmcVector<double> vc_v;

	Reconstruct(vc_v, vc_a);
	
	int i,j;

	for (i=0; i<vc_v.size(); i++)
		vc_v(i) = vc_v(i)*(dmax-dmin);

	if (ras_image.rows()==HeightOfPCA() && ras_image.cols()==WidthOfPCA())
	{
		for (j=0; j<vc_v.size(); j++)
			ras_image(j) = vc_v(j);
	}
	else
	{
		Raster<double>	ras_reconstruct;

		ras_reconstruct.SetDimension(HeightOfPCA(), WidthOfPCA());
		
		for (j=0; j<vc_v.size(); j++)
			ras_reconstruct(j) = vc_v(j);

		ras_reconstruct.ResizeTo(ras_image, INTRP_BILINEAR);
	}
}


void MCMC_FacePCA::SubtractMean(mcmcVector<double> &vc_v, mcmcVector<double> &vc_x)
{
	assert(vc_x.size() == SizeOfPCA());

	vc_v.SetColVector(SizeOfPCA());
	for (int i=0; i<SizeOfPCA(); i++)
		vc_v(i) = vc_x(i)-vc_mean(i);
}



double MCMC_FacePCA::Prob(mcmcVector<double> &vc_x, mcmcVector<double> *pvc_a)
{
	// subtract the mean
	mcmcVector<double> vc_v;

	SubtractMean(vc_v, vc_x);

	double d1 = pow(vc_v.Norm(2), 2.0);
		
	mcmcVector<double> vc_a1;

	// projection to the PCA
	Projection(vc_a1, vc_x);

	int		i;
	double	d2,d3,sum,val;
	
	// initialization
	d2 = 0.0;
	d3 = 0.0;
	sum = 0.0;

	for (i = 0; i < NumComponents(); i++)
	{
		val = de / (vc_sigmas(i) + de);
		d3 += log(val);
		val = vc_sigmas(i) / (vc_sigmas(i) + de);
		d2 += (vc_a1(i) * vc_a1(i) * val);
		sum += vc_sigmas(i);
	}


	d3/=2;
	double d4 = -SizeOfPCA() * log(PI * 2 * de) / 2;

	double p = (d4 + d3 + d2/2/de - d1/2/de);
	
	if (pvc_a != NULL)
		*pvc_a = vc_a1;

	return p;
}

//
//	Get best face probability in the region and synthesize
//
double MCMC_FacePCA::MatchFaceModel(Raster<float> &image_input, irmm_Rect rect_to_test, MCMC_FaceRect &face)
{
	int				i,j,k,size_test;
	Raster<float>	ras_face, ras_normalized;
	double			dtemp,e;
	irmm_Rect		rect_try;
	mcmcVector<double>	vc_v,vc_a_temp;
	float			intensity_min, intensity_max;
	int			m;
	
	// initialization
	if (rect_to_test.Height() < rect_to_test.Width())
		size_test = rect_to_test.Height();
	else
		size_test = rect_to_test.Width();

	face.energy = 1e10;

	// set sizes
	vc_v.SetColVector(SizeOfPCA());
	ras_normalized.SetDimension(HeightOfPCA(), WidthOfPCA());
	face.rect_face = rect_to_test;

	for (k=0; k<2; k++)
	{
		if (size_test < 10) break;

		for (j = rect_to_test.top; j <= rect_to_test.bottom+1 - size_test; j++)
			for (i = rect_to_test.left; i<= rect_to_test.right+1 - size_test; i++)
			{
				rect_try.Set(i, j, i+size_test-1, j+size_test-1);

				ras_face.SetDimension(rect_try.Height(), rect_try.Width());

				image_input.CopySubImage(rect_try, ras_face);

				ras_face.ResizeTo(ras_normalized, INTRP_BILINEAR);
				
				// normalize pixel values
				ras_normalized.GetMinMax(intensity_min, intensity_max);
				dtemp = 1.0/(intensity_max-intensity_min);

				assert(vc_v.size() == ras_normalized.size());
				for (m=0; m<ras_normalized.size(); m++)
					vc_v(m) = dtemp*ras_normalized.Data(m);

				e = 2000 - Prob(vc_v, &vc_a_temp);
				if (e < 0) e = 0;
				
				if (e < face.energy)
				{
					
					face.energy = e;

					face.rect_face = rect_try;
					face.vc_a_face = vc_a_temp;

					face.dmin = intensity_min;
					face.dmax = intensity_max;
				}
			}
		
		size_test = IntClose(1.0/1.1*size_test);
	}

	return face.energy;
}

////////////////////////////////////////
////////////////////////////////////////

// construction & destruction
MCMC_FaceParticles::MCMC_FaceParticles()
{
	face_pca.LoadFromFile("face.pca");
}

MCMC_FaceParticles::~MCMC_FaceParticles()
{
}

void MCMC_FaceParticles::Empty(void)
{
	list_center_rect.Empty();

	vc_biggest_center_rect.SetDimension(0,0);
	vc_precomputed_facemodels.SetDimension(0,0);
}

// operations

void MCMC_FaceParticles::PrecomputeFaceModels(CFaceDetector &face_detector, Raster<float> &ras_intensity)
{
	Empty();
	
	if (face_detector.NumDetectedRect() <= 0)
		return;

	// get all the rects
	bool			bgo,bnear;
	RectOfInterest	roi;
	int				i,j,k;
	McPoint2D		pt_center_existing,pt_center_current;
	
//	Raster<double>	ras_syn;
//	Raster<uchar>	ras_display;
//	char str[200];

	vc_precomputed_facemodels.SetDimension(1, face_detector.NumDetectedRect());
	k = 0;
	bgo = face_detector.GetFirstDetectedRect(roi);
	while (bgo)
	{
		// precompute the face models
		face_pca.MatchFaceModel(ras_intensity, roi.Rect(), vc_precomputed_facemodels(k));
		
/*		ras_syn.SetDimension(vc_precomputed_facemodels(k).rect_face.Height(), vc_precomputed_facemodels(k).rect_face.Width());
		face_pca.Reconstruct(ras_syn, vc_precomputed_facemodels(k));
		sprintf(str, "%d.tif", k);
		CopyRaster(ras_display, ras_syn);
		SaveRaster2ImageFile(ras_display, str);
*/
		vc_precomputed_facemodels(k).rect_face.Center(j, i);
		pt_center_current.Set(j, i);
		
		// put it into the center list
		bnear = false;
		bgo = list_center_rect.GetFirst(pt_center_existing);
		while (bgo)
		{
			if (Distance(pt_center_current, pt_center_existing) < 15.0)
			{
				bnear = true;
				break;
			}
			bgo = list_center_rect.GetNext(pt_center_existing);
		}
		if (!bnear)
			list_center_rect.Insert(pt_center_current);


		bgo = face_detector.GetNextDetectedRect(roi);
		k++;
	}

	// the biggest rect centered at each proposal position
	mcmcVector<int>	vc_size_temp;
	int			size_temp;

	vc_biggest_center_rect.SetDimension(1, list_center_rect.Num());
	vc_size_temp.SetRowVector(vc_biggest_center_rect.size());
	vc_size_temp.InitValue(0);
	for (k=0; k<NumFaceRect(); k++)
	{
		vc_precomputed_facemodels(k).rect_face.Center(j, i);
		pt_center_current.Set(j, i);

		j = 0;
		bgo = list_center_rect.GetFirst(pt_center_existing);
		while (bgo)
		{
			if (Distance(pt_center_current, pt_center_existing) < 15.0)
			{
				size_temp = vc_precomputed_facemodels(k).rect_face.Area();
				if (size_temp > vc_size_temp(j))
				{
					vc_size_temp(j) = size_temp;
					vc_biggest_center_rect(j) = vc_precomputed_facemodels(k).rect_face;
				}
			}
			j++;
			bgo = list_center_rect.GetNext(pt_center_existing);
		}
	}
}

int	MCMC_FaceParticles::BestFitRect(double darea_region, irmm_Rect bound_region)
{
	McPoint2D	pt_center,pt_center_rect;
	int			i_rect;
	bool		bnear,bgo;

	i_rect = -1;
	bnear = false;
	
	double dtemp=IntToDouble(bound_region.Area());
	if (darea_region/dtemp < 0.40)
		return -1;	// not a compact region

	pt_center.x = (bound_region.top+bound_region.bottom)/2.0;
	pt_center.y = (bound_region.left+bound_region.right)/2.0;

	bgo = list_center_rect.GetFirst(pt_center_rect);
	while (bgo)
	{
		if (Distance(pt_center_rect, pt_center) < 20.0)
		{
			bnear = true;
			break;
		}
		bgo = list_center_rect.GetNext(pt_center_rect);
	}

	if (!bnear)
		return -1;	// too far from the face center from bottom-up processes

	int	larea_biggest=-1,loverlapping;

	for (int i=0; i<NumFaceRect(); i++)
	{
		loverlapping = bound_region.OverlappingArea(vc_precomputed_facemodels(i).rect_face);
		double aa=IntToDouble(loverlapping)/IntToDouble(vc_precomputed_facemodels(i).rect_face.Area());
		if (aa > 0.65 &&
			IntToDouble(loverlapping)/darea_region > 0.40)
		{
			if (loverlapping > larea_biggest)
			{
				i_rect = i;
				larea_biggest = loverlapping;
			}
		}
	}
	return i_rect;
}

bool MCMC_FaceParticles::InFaceRects(const int row, const int col)
{
	for (int i=0; i<vc_biggest_center_rect.size(); i++)
	{
		if (vc_biggest_center_rect(i).InRect(row, col))
			return true;
	}
	return false;
}

void MCMC_FaceParticles::ComputeFaceModel(LinkList<PMCMC_PixelSet> &list_domain, Raster<float> &ras_intensity, MCMC_FaceRect &face_rect)
{
	double		darea_region;
	irmm_Rect	bound;
	
	irmm_Rect	bound_temp;
	bool		bgo_domain;
	MCMC_PixelSet	*ppixel_set;

	// to compute the bound of the region
	bound.Set(9999,9999,-9999,-9999);
	darea_region = 0.0;

	bgo_domain = list_domain.GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bound_temp = ppixel_set->Bound();
			darea_region += ppixel_set->Num();

			bound.SetBound(bound_temp.left, bound_temp.top);
			bound.SetBound(bound_temp.right, bound_temp.bottom);
		}

		bgo_domain = list_domain.GetNext(ppixel_set);
	}

	// fit existing rectangle?
	int	irect;

	irect = BestFitRect(darea_region, bound);
	if (irect >= 0)
		face_rect = vc_precomputed_facemodels(irect);
	else
		face_pca.MatchFaceModel(ras_intensity, bound, face_rect);
}



////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
//        Models for face

MCMC_FaceParticles *pface_saliency_map;


//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
MCMC_Intensity_Face::MCMC_Intensity_Face()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_FACE;

	type = MODEL_INTENSITY_FACE;
}

MCMC_Intensity_Face::~MCMC_Intensity_Face()
{
	Free();
}

void MCMC_Intensity_Face::Free(void)
{
}

// operators
MCMC_Intensity_Face & MCMC_Intensity_Face::operator=(const MCMC_Intensity_Face &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);
		
		surface			= x.surface;
		intensity_mean	= x.intensity_mean;
		intensity_var	= x.intensity_var;
		bin_energy		= x.bin_energy;

		face_rect		= x.face_rect;

		pface_particles	= x.pface_particles;
	}

	return *this;
}

MCMC_Intensity_Face & MCMC_Intensity_Face::operator=(MCMC_Model *pa)
{
	if (pa!=NULL && this!=pa && pa->Type()==Type())
		*((MCMC_Intensity_Face*) this) = *((MCMC_Intensity_Face*) pa);

	return *this;
}

int MCMC_Intensity_Face::TotalCodingBits(void)
{
	int bits;

	bits =	surface.TotalCodingBits();		// control points
	bits += sizeof(float)*39;

	return bits;
}


void MCMC_Intensity_Face::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);
	
	strm.write((char*) &intensity_mean, sizeof(intensity_mean));
	strm.write((char*)&intensity_var, sizeof(intensity_var));
}

void MCMC_Intensity_Face::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);

	strm.read((char*)&intensity_mean, sizeof(intensity_mean));
	strm.read((char*)&intensity_var, sizeof(intensity_var));
}

MCMC_Model* MCMC_Intensity_Face::CopyTo(void)
{
	MCMC_Intensity_Face *preturn;

	preturn = new MCMC_Intensity_Face();
	(*preturn) = (*this);
	
	return ((MCMC_Model*)preturn);
}

bool MCMC_Intensity_Face::MeanFromPixelSet(void)
{
	// compute the intensity means
	int				i,j;
	mcmcVector<double>	vc_temp(3);
	mcmcMatrix<int>		mx_temp;

	assert(pface_particles != NULL && plist_domain != NULL);

	pface_particles->ComputeFaceModel(*plist_domain, prespFilter->ras_intensity, face_rect);

	// compute the surface first
	surface.Initialize(*plist_domain, pmx_taken, prespFilter->ras_intensity);
	surface.Estimate();

	intensity_mean = 0.0;
	nsize = 0;
	
	MCMC_PixelSet	*ppixel_set;
	bool			bgo_domain,bgo;

	// compute the major direction
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (!face_rect.InRect(j,i))
				{
					// get intensity
					intensity_mean += prespFilter->ImageIntensity(j,i);
					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize > 0)
		intensity_mean = intensity_mean/nsize;

	return true;
}

bool MCMC_Intensity_Face::VarianceFromPixelSet(void)
{
	int		i,j;
	double	din,num_small,small_var,num_big,num_count,intensity_onsurface;
	int	num_pixels_nonfaces;
	
	num_pixels_nonfaces = 0;

	// initialization
	nsize = 0;
	intensity_var = 0.0;
	num_small = 0.0;
	small_var = 0.0;
	num_big = 0.0;
	num_count = 0.0;
	
	assert(plist_domain != NULL);

	MCMC_PixelSet	*ppixel_set;
	bool			bgo_domain,bgo;

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (!face_rect.InRect(j,i))
				{
					// projection
					din = prespFilter->ImageIntensity(j,i);
					intensity_onsurface = surface.SurfaceValue(j,i);

					// for robust variance
					if (fabs(din-intensity_onsurface)<5.0)
						num_small += 1.0;
					if (fabs(din-intensity_onsurface)>35.0)
						num_big += 1.0;
					else
					{
						small_var += pow(din-intensity_onsurface,2.0);
						num_count += 1.0;
					}

					// variance of intensity
					intensity_var += pow(din-intensity_onsurface,2.0);
					num_pixels_nonfaces++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
			nsize++;
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (num_pixels_nonfaces>1)
	{
		if (num_pixels_nonfaces>1000 && num_small/((double)num_pixels_nonfaces)>0.90 && num_big/((double)num_pixels_nonfaces)<0.04)
			intensity_var = small_var/num_count;
		else
			intensity_var = intensity_var/(num_pixels_nonfaces-1);
	}
	else
		intensity_var = 1.0;

	if (intensity_var < 5.0)
		intensity_var = 5.0;
		
	// compute quantized energy
	Histogram	hist_temp;
	double		value=0.0,prob,dis;

	bin_energy.Set(0.0,INTENSITY_HIGH,INTENSITY_DELTA);
	hist_temp.Set(0.0,INTENSITY_HIGH,INTENSITY_DELTA);

	while (value<INTENSITY_HIGH+INTENSITY_DELTA)
	{
		dis = pow(value,2.0)/intensity_var/2.0;
		prob = ProbFromEnergy(RobustLoglik(dis)+log(sqrt(2.0*PI))+log(intensity_var)*0.5);
		hist_temp.AddBin(value,prob);
		value += 1.0;
	}
	
	// hist_temp should not be normlized
	for (i=0; i<hist_temp.count; i++)
	{
		if (hist_temp.pbins[i]>0.0)
			bin_energy.pbins[i] = -log(hist_temp.pbins[i]);
		else
			bin_energy.pbins[i] = 60;
	}

	return true;
}

bool MCMC_Intensity_Face::LogLikelihoodFromPixelSet(void)
{
	int				i,j;
	double			dis_uc,dis,din,dshift;
	double			darea_rect,num_noneface_pixels;

	darea_rect = (double)face_rect.rect_face.Area();

	// initialization
	nsize	= 0;
	log_pos_prob = 0.0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;
	num_noneface_pixels = 0;

	double	real_size;
	real_size = 0.0;

	assert(plist_domain != NULL);

	MCMC_PixelSet	*ppixel_set;
	bool			bgo_domain,bgo;

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				real_size += 1.0;
				
				if (!face_rect.InRect(j,i))
				{
					// LogLikelihood in intensity
					din = prespFilter->ImageRegionLayerIntensity(j,i);
					dshift = fabs(din-surface.SurfaceValue(j,i));

					dis = pow(dshift,2.0)/intensity_var/2.0;
					log_pos_prob += bin_energy.GetBin(dshift);
						
					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(bin_energy.GetBin(dshift)+Complexity());
					crying_voice += dis_uc;

					num_noneface_pixels += 1.0;
				}
				nsize++;
	
				bgo = ppixel_set->GetNext(j, i);
			}
		}

		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;
		
	log_pos_prob += pow(num_noneface_pixels,COMPLEXITY_POW)*scale_factor+
					   pow((double)nsize-num_noneface_pixels,COMPLEXITY_POW)*scale_factor;
	log_pos_prob += num_noneface_pixels*1.0;
	log_pos_prob += ((double)nsize)*Complexity();
	log_pos_prob += face_rect.energy*1.0/1.5*(nsize-num_noneface_pixels)/darea_rect;

	crying_voice = pow(crying_voice, 0.6);

	return true;
}

void MCMC_Intensity_Face::UpdateLogLikelihood(const int nlabel, MCMC_PixelSet &list_pixels_old, MCMC_PixelSet &list_pixels_new)
{
	//register	i,j;
	int i, j;
	double		dshift,din,dtemp,dcomplexity_temp;
	bool		bgo;

	if (list_pixels_old.Num()<=0 && list_pixels_new.Num()<=0)
		return;
	
	log_pos_prob -= pow((double)nsize,COMPLEXITY_POW)*scale_factor;

	dcomplexity_temp = Complexity();
	// delete the old pixels
	bgo = list_pixels_old.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dshift = fabs(din-surface.SurfaceValue(j,i));
			dtemp = bin_energy.GetBin(dshift);

			log_pos_prob += dtemp+dcomplexity_temp;
			log_lik += dtemp;

			nsize++;
		}
		bgo = list_pixels_old.GetNext(j, i);
	}

	bgo = list_pixels_new.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dshift = fabs(din-surface.SurfaceValue(j,i));
			dtemp = bin_energy.GetBin(dshift);

			log_pos_prob = log_pos_prob-(dtemp+dcomplexity_temp);
			log_lik -= dtemp;

			nsize--;
		}
		bgo = list_pixels_new.GetNext(j, i);
	}

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
}


void MCMC_Intensity_Face::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->ImageIntensity(row_pixel,col_pixel);
}

double MCMC_Intensity_Face::GetDirectDis(mcmcVector<double> &vc_intensity)
{
	double	din;
	
	din = vc_intensity(0);

	return (intensity_mean-din)*(intensity_mean-din); // give lower prevelige over intensity
}


double MCMC_Intensity_Face::GetPdf(void)
{
	return 1.0;
}

void MCMC_Intensity_Face::PrepreForSynthesize(void)
{
	rf_i.Set(0.0,sqrt(intensity_var),1);

	assert(pface_particles != NULL);

	ras_reconstruct.SetDimension(face_rect.rect_face.Height(), face_rect.rect_face.Width());
	pface_particles->Reconstruct(ras_reconstruct, face_rect);
}

bool MCMC_Intensity_Face::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a.SetRowVector(3);

	if (!face_rect.InRect(row, col))
		vc_a(0) = surface.SurfaceValue(row,col)+rf_i.Sampling()/1.5;
	else
	{
		int j,i;

		j = row-face_rect.rect_face.top;
		i = col-face_rect.rect_face.left;
		vc_a(0) = ras_reconstruct(j, i)+rf_i.Sampling()/5;
	}

	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	for (int i=0; i<vc_a.size(); i++)
	{
		if (vc_a(i)<0) vc_a(i) = 0;
		if (vc_a(i)>255) vc_a(i) = 255;
	}

	return true; // in RGB space
}


bool MCMC_Intensity_Face::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a.SetRowVector(3);

	if (!face_rect.InRect(row, col))
		vc_a(0) = surface.SurfaceValue(row,col);
	else
	{
		int j,i;

		j = row-face_rect.rect_face.top;
		i = col-face_rect.rect_face.left;
		vc_a(0) = ras_reconstruct(j, i);
	}

	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	for (int i=0; i<vc_a.size(); i++)
	{
		if (vc_a(i)<0) vc_a(i) = 0;
		if (vc_a(i)>255) vc_a(i) = 255;
	}

	return true; // in RGB space
}


// get the log of the posterior probability at (row,col)
double MCMC_Intensity_Face::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	assert(prespFilter!=NULL);
	
	return 0.0;
/*
	int				i,j;
	double			dis=0.0,count=0.0,dshift,darea_rect;
	MCMC_FaceRect	*pface_rect;

	pface_rect = pface_saliency_map->GetFaceRect(i_rect);
	if (pface_rect != NULL)
		darea_rect = (double)pface_rect->rect_bound.Area();

	for (j=row-1; j<=row; j++)
		for (i=col-1; i<=col; i++)
			if (prespFilter->Valid(j,i) && (pmx_occupied==NULL || (pmx_occupied!=NULL && pmx_occupied->Data(j,i)==0)))
			{
				if (pface_rect != NULL && pface_rect->rect_bound.InRect(j, i))
					dis += pface_rect->energy/1.5/darea_rect;
				else
				{
					dshift = fabs(prespFilter->ImageRegionLayerIntensity(j,i)-surface.SurfaceValue(row,col));

					dis += bin_energy.GetBin(dshift);
				}
				count = count+1.0;
			}

	if (count < 2.0 && pmx_occupied != NULL)
	{
		dis = 0.0;
		count = 0.0;

		for (j=row-1; j<=row+1; j++)
			for (i=col-1; i<=col+1; i++)
				if (prespFilter->Valid(j,i) && pmx_occupied->Data(j,i)==0)
				{
					dshift = fabs(prespFilter->ImageRegionLayerIntensity(j,i)-surface.SurfaceValue(row,col));

					dis += bin_energy.GetBin(dshift);
					count = count+1.0;
				}
	}

	if (count>0)
		dis = dis/count;

	return dis+Complexity();
*/
}


void MCMC_Intensity_Face::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "%.2f(%.5f) surface %.1f:%.1f", log_pos_prob,log_lik/nsize,intensity_mean,intensity_var);
}
